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ABSTRACT 



o 

^ Context. Modelling the gaseous component of the interstellar medium (ISM) by Smoothed Particles Hydrodynamics in N-Body 
C/3 , simulations (NB-TSPH) is still very crude when compared to the complex real situation. In the real ISM, many different and 
> almost physically decoupled components (phases) coexist for long periods of time, and since they spread over wide ranges of 
' density and temperature, they cannot be correctly represented by a unique continuous fluid. This would certainly influence star 
, formation which is thought to take place in clumps of cold, dense, molecular clouds, embedded in a warmer, neutral medium, 
^ ' that are almost freely moving throughout the tenuous hot ISM. Therefore, assuming that star formation is simply related to 
OO ' the gas content without specifying the component in which this is both observed and expected to occur, may not be physically 
. sounded. 

Aims. In this paper we seek to improve upon this aspect by considering a multi-phase representation of the ISM in NB-TSPH 
simulations of galaxy formation and evolution with particular attention to the case of early-type galaxies. 

Methods. Cold gas clouds are described by the so-called sticky particles algorithm. They can freely move throughout the hot 
' ISM medium; stars form within these clouds and the mass exchange among the three baryonic phases (hot gas, cold clouds, 
t , stars) is governed by radiative and Compton cooling and energy feedback by supernova (SN) explosions, stellar winds, and UV 
' radiation. We also consider thermal conduction, cloud-cloud collisions, and chemical enrichment. 
^ ' Results. Our model agrees with and improves upon previous studies on the same subject. The results for the star formation rate 
• i-H , are very promising and agree with recent observational data on early-type galaxies. 

Conclusions. These models lend further support to the revised monolithic scheme of galaxy formation, which has recently been 
also strengthened by high redshift data leading to the so-called downsizing and top-down scenarios. 

. - _ t 
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1. Introduction 

Since almost two decades, N-Body (NB) simulations are 
the best tool we have to study complex systems such 
as galaxies and clusters. A number of different recipes 
have been proposed to model the gaseous component of 
such systems. The best known is probably the Smoothe d 
Particle Hydrodynamics (SPH) (|Lucvlll977t lBenj|l990l ). 
In brief, it describes the gaseous fluid with a finite num- 
ber of particles, which are followed in their Lagrangian 
evolution, and evaluates the values of the physical quan- 
tities in each point of the real space by smoothing each 
particle through a kernel function, peaked in its position. 
Together with the equations of hydrodynamics it allows 
us to study mechanical heating and shocks; coupled to 



gravity described by the well known Tree code (T) by 
Barnes fc (|l986f) . the NB-TSPH method has proved 
to be very efficient. 

Nevertheless, it suffers some limitations. In particu- 
lar, it seems not to be accurate enough to follow in de- 
tail the extreme complexity of the structure and physical 
processes occurring in a galaxy, which include radiative 
cooling, star forma tion, feedback from Supernovae, and 
chemical evolution (lKatzlll992l: iMihos fc Hernauistl ll994 



Chiosi fc Garrard 120021; iMerhn fc Chiosil |2006| ). These 
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processes, apart from requiring sufficient computational 
power in order to achieve the necessary resolution, also 
need some kind of different, ad hoc descriptions for the 
various components. The subject is still at infancy, which 
in part explains the many uncertainties currently affect- 
ing the theoretical models. Star formation is usually ex- 
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pressed by the empiric law of ISchmidtl ([l95^, in which 
the rate is proportional to the gas density, often cou- 
pled with so me conditions rem iniscen t of the J eans sta- 
bility (see e.g. iGerritsen fc Icke 1997; Carraro et al. 1998; 



Im arri fc White! [iooi . How to share chemical elements 
ejected by a star particle with the surrounding medium 
is another point of unc ertainty. Also in this case the sta- 



Buonomo et al. 200(J) . If the local gas density obeys the 



above conditions, gas particles are supposed to be turned 
into star particles. In the so called sticky-particles al- 
gorit hm ( Levinson fc Roberts! 1981 : Noguchi fc Ishibashi 
19861) star formation is considered to be enhanced by 
cloud-cloud collisions. 

Different recipes have been tested to introduce the 
complicated physics of star formation into numerical 
codes. The most physically grounded, although expen- 
sive in computational load, simply consists in creating 
new star particles of different mass according to a cho- 
sen initial mass function (IMF). Obviously this causes a 
huge growth of the total number of particles, which is well 
beyond present-day technical capability. The situation is 
more viable if, assumed a certain mass resolution (initial 
total baryonic mass in form of gas divided by the number 
of mass particles used in the simulation), a single star par- 
ticle is conceived as an assembly of real stars whose mass 
distribution follows a given IMF. In this context, the total 
number of baryonic particles does not diverge. Therefore 
a single gas particle can be turned into a star particle 
with no serious consequences for the NB-TSPH ecology. 
What happens inside a star particle is relevant only for 
the energy and chemical yields but not for the NB-TSPH 
description as such. This scheme has been follow e d in a 
number of studies. For instance, lJungwiert et al.! (|200ll ) 



considers the gradual formation of stars inside gas par- 
ticles, which eventually become star particles when the 
hidden star mass exceeds a certain threshold value. Along 
this line of th ought is the probabilis ti c descript i on of star 
formation by iBuonomo et al. ( 2000l ) ; Lia et al. ( 2002! ) , in 
which gas particles are turned into star particles accord- 
ing to the Schmidt law however interpreted in statistical 
sense, so that at any time there is a certain probability 
that a gas particle becomes a star particle. This way of 
looking at the star formation process is particularly ad- 
vantageous in numerical simulations as it keeps constant 
the total number of particles, thus saving lots of compu- 
tational time. 

Owing to the current mass resolution of numerical sim- 
ulations, the star particles actually have the mass size of 
a star cluster, in which the real stars are born in a short 
burst of intense local activity with homogeneous age and 
chemical composition. Therefore each star particle can 
be approximated to a single stellar population (SSP), in 
which the stars distribute in mass according to some IMF. 

The influence of the IMF describing the stellar con- 
tent of each star particle becomes soon evident when the 
energy and chemical feedback by supernova (SN) explo- 
sions and/or stellar winds to the surrounding medium are 
taken into account. Several unsettled questions come soon 
into play. Among others, it is still not clear whether the 
SN energy feedback should be only treated as "thermal" 
or if some kinetic effects have to be considered (see e.g. 



2002!) turns out to be 



tistical description by iLia et al" 
particularly convenient. 

In current NB-TSPH models, the gas component is 
considered as an unique fluid regardless of the observa- 
tional fact that in the gas content of real galaxies several 
different phases are present: hot and cold, tenuous and 
dense, ionised, neutral and molecular gas. Furthermore, 
in the same models star formation and its efficiency are 
related to the total gas content with no attention to the ob- 
servational fact that stars are formed only in dense molec- 
ular clouds. In this paper, we present a multi-phase de- 
scription of the ISM (hot gas and cold clouds with molec- 
ular cores inside) in which star formation can occur only 
in the cold component thus bringing the simulations closer 
to reality. The new prescription for the ISM is imple- 
mented in G alDyn, the Pa dua NB-TSPH code originally 
developed by Carraro et al.! d 19981) and later improved by 
iLia et al.! (|2002l ) and lMerlin fc Chiosi! (|20n6h . 

The layout of the paper is as follows. In Sect. [2] we 
present in some details the physical problem and briefly 
summarize our prescription for the ISM and how it is im- 
plemented in the numerical code. In Sect.[3]we test the new 
method simulating the formation of an early-type galaxy 
and present the results. In Sect. H] we discuss the results, 
compare them to some recent observational data on galaxy 
formation and evolution, and draw some general conclu- 
sions. 



2. Multi-phase ISM 

2.1. The physical context 

Since the early sixties, the ISM in the Milky Way is known 
not to be uniform on small scales, but to contain strong lo- 
cal density enhancements, so that very dense and cold re- 
gions, in form of atomic and molecular " clouds " , appear to 



be em bedded in a diffuse hot, ionised medium ([Field et al 
19691 ). Some years later, with the increasing understand- 



ing of the importance of SN explosions, another compo- 
nent was introduced, i.e. the hot and very tenuous corona l 
gas expelled by exploding stars i McKee fc Ostrikeilll977 ). 

Therefore, it is nowadays widely accepted that the ISM 
of the Milky Way is made of several components with dif- 
ferent temperatures and densities. At least four, almost 
decoupled, gaseous phases are present: (i) the hot, highly 
or fully ionised, tenuous medium, with temperatures of 
some tens of thousands of Kelvin degrees and number 
densities of about 0.1 cm~'^; (ii) the cold, dense, clumpy 
molecular medium, with temperatures of only a few Kelvin 
degrees and number densities from thousands up to mil- 
lions of atoms per cm^; (iii) a cloud-embedding, warm, 
neutral phase, with temperatures from a few thousands 
down to a few hundreds Kelvin degrees, and number den- 
sity of a few atoms per cm^; (iv) and an extremely hot and 
rarefied coronal gas expelled by the SN explosions, with 
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temperatures of millions of Kelvin degrees and very low 
number densities (some 0.01 atoms per cm'^). Of course, 
this scheme is still a simple one, because the ISM contains 
other local features (HII regions, for instance, which are 
regions of ionised and hot but very dense material) . 

In the Milky Way, these different phases fill the in- 
terstellar space in a way that is inversely proportional to 
their density. Most of the space is filled by the coronal and 
hot gas, whereas only a minor fraction of it is filled by the 
cold clouds. The cold clouds, on the other hand, contain 
the largest fractional mass of the ISM. Owing to such a 
strong density contrast, the cold clouds are not strongly 
affected by drag interactio ns with the tenuous phase (see 



Levinson fc Robert^ 1981 . , and below). 

Given this situation, it seems clear that the SPH algo- 
rithm, which is fully adequate to describe a single-phase 
fluid even in situations of shock, cannot be considered the 
best tool to study such a complex mixture of different 
phases in mechanical and thermodynamical equilibrium. 
A well known problem of the SPH method, for instance, 
is the so called overcooling. Considering the whole mass of 
gas as a single fluid, a particle can have neighbouring par- 
ticles with much higher densities. These should actually 
be considered as part of another decoupled phase. In the 
SPH formalism they actually end up to dominate the local 
density of the particle under consideration, thus leading to 
overestimate the density and con sequently to a strong er, 
un-physical radiative cooling (see Marri &: White! 20031 ). 

Furthermore, the correct equation of state relating 
pressure, density, temperature, and chemical composition 
of the gas should properly include the molecular compo- 
nent. In most numerical NB-TSPH codes including our 
own, the equation of state is that for an ideal mono-atomic 
gas, which obviously does not adequately describe the 
molecular component. Finally, self-gravitating molecular 
clouds are kn own to obey particular empirical laws (e.g. 
Larson 198l[ ). Likely, they should be taken into account 
in the grand-design of star formation. 

Over the years, many attempts have been made 
to solve this problem and to find a physically 
sounded recipe for s tar fo rmation in molecular clouds. 
Levinson fc Robert^ ( 1981 ) developed the so called sticky 
particles algorithm, in which the ISM is modelled only by 
means of clouds freely moving across t h e surr ounding void 
space. More recently. iMarri fc White ( 2003f ) suggested a 
simple modification of the neighbour searching algorithm 
in SPH: particles of very different density and temperature 
are discarded in the neighbour sam pling. Other strate- 
gies h a ve been recently pr opos ed by Semelin fc Combe"^ 
(|2002[) . iHarfst et all (|2006') and B ooth et al.l (|2007t ). trv- 
ing to simulate the evolution of disk galaxies. 

A description of the ISM as detailed as possible is re- 
quired also to study galaxy formation in cosmological con- 
text. Even if there are no evidences that the primordial 
ISM was similar to the present-day ISM in the Milky Way, 
nor the contrary, it is likely that the ISM always had an 
important role. In particular, the duration and intensity 
of the star formation activity should be strongly depen- 



dent on how different components of matter interact each 
other. Along this line of thou ght are the attempts made by 
Springel fc Hernguistl (|2003l ) to reproduce the complexity 
of the ISM in cosmological context. 

2.2. The multi-phase prescription 

Starting from these considerations and the above men- 
tioned studies, we tried to elaborate a recipe for the var- 
ious components of the ISM and their relationships. In 
brief, hot gas is slowly turned into cold clouds and star 
formation occurs in these latter. Viceversa the energy de- 
posit by dying stars heats up the cold gas. The key points 
of our recipe are the listed below. 

(1) We assume that the SPH method can reasonably 
describe the hot phase of the ISM, considering as "hot" 
the gas particles which are hotter and less dense than suit- 
able temperature and density thresholds (i.e. both the hot 
ionised gas and the coronal gas). 

(2) The cold phase (molecular clouds), in which star 
formation occurs, is described by a different method, based 
on the fact that individual particles can roughly represent 
single clumps of cold material. A giant molecular cloud 
can reach a mass of some 10^ M0. Considering also the 
embedding neutral or slightly ionised warm material, the 
total mass of a cold particle can be higher than this value. 
As discussed below, the mass resolution of our fiducial 
model is close to 6 x 10^ Mq thus making our description 
quite realistic. 

(3) Following levinson fc Roberts! (|l98l[ ). we assume 
that the cold clumps can move almost freely throughout 
the hot phase. Pressure drag and friction are certainly 
present, but we consider them negligible for our pur poseqj. 

(4) As suggested bv ISemelin fc CombesI ( 20021 ). a hot 
gas particle does not instantaneously turn into a cold cloud 
particle when its temperature falls below the threshold 
value because letting this to occur would originate too 
sharp a transition. Therefore, we prefer to use a smooth 
procedure named gathering method. Each gas particle is 
thought to be made of gas at different temperatures dis- 
tributed according to a gaussian law, truncated at 2(7 from 
the central peak determined by the SPH value. To some 
extent, this mirrors the situation of the SPH technique, 
in which each particle has properties that result from av- 
eraging those of nearby particles, so that only groups of 
particles can be considered as representative of the dis- 
tribution of matter in a given region. The (Tco/d of this 
quasi-gaussian distribution is a free parameter. Each par- 
ticle is then assigned a certain fraction of cold gas, i.e. the 
amount of gas in the quasi-gaussian distribution whose 
temperature falls below the threshold value Tcoid- At each 
time step, a particle tries to "acquire" some cold gas from 
its SPH neighbours until a single cold particle of mass 
equal to that of the mother particle is formed. When this 
is achieved, a new " cold" particle is created from the " hot" 



' See e.g. IShu et all (| 19721 ) for an estimate of ram pressure 
due to the ambient gas. 
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one. The cold particle is separated from the surrounding 
hot medium, provided it belongs to a convergent flow (i.e. 
its velocity divergence is negative), and its cooling time is 
shorter than its dynamical time, whose evaluation take s 
Dark Matter around into account ( Buonomo et al.|[200Cll ). 
The temperature Thi of the new cold particle is deter- 
mined by the condition that during the creation process 
the total thermal energy of the group of particles remains 
constant, together with the requirement that neighbours 
have no longer cold gas (their temperature distributions 
accordingly vary). If the temperature of a gaseous hot par- 
ticle is limited by the Jeans threshold (see Sect. 13.11 point 
3), its acoid is changed so that it can provide some small 
fraction of cold gas anyway. The limitation imposed by 
the Jeans threshold is indeed a numerical artifact due to 
resolution limits, and there are no real physical reasons to 
exclude the presence of a small fraction of cold gas within 
the total gaseous mass of the particle. The acoid is ad- 
justed in such a way that the lower limit of the truncated 
gaussian distribution corresponds to zero temperature. 

Finally, we assign the new cold particle the initial den- 
sity obtained by interpolating with a cubic polynomial 
fit the rel ation between temperature and number density 
given by IScheffler fc Elsaessei ( 1987 ) for cold clouds in 
pressure equilibrium with the surrounding medium: 



log[nHi] = -0M3x^ + S.ix"^ - 24.66a; + 23.5, (1) 
where x = log[THi]- The relation in shown in Fig. [1] 



causing fragmentation in cold, extremely dense cores of 
molecular gas. N umerical simulations of the impact be- 
tween two clouds ( Marinho Sz Lepindl200d iMarinho et al 



log[T] 

Fig. 1. Number density versus temperature of atomic 
Hydrogen in pressure equilibrium. 

(5) In addition to gravity, the only interactions that 
cold particles can feel are the mutual inelastic collisions 
(jLevinson fc Roberts! [l98ll ). 

Collisio ns between clouds are known to trigger star 
formation ( Loren 119761 : iHasegawa et al.|[l994[) . due to the 
creation of a hot shock front which rapidly cools down 



[2001.) showed how the results of such impacts depend on 
many parameters, often acting well below our resolution 
scale. Anyway, based on those simulations as well as on 
the cited observational evidences, we can envision the fol- 
lowing scenario. First of all, two clouds are considered to 
be colliding if their mutual distance is shorter then the 
sum of the radii of each particle. A rigourous description 
of cloud-cloud collisions would require the use of the im- 
pact parameter instead of the distance. However, in the 
framework of the statistical approach used for star forma- 
tion and gas restitution (see below) using the distance is 
adequate to our purposes. The radius of a particle is de- 
rived from its mass to total density ratio and the assump- 
tion of spherical symmetry. The total density and mass 
include both the neutral and the molecular components. 
Second, the densities of the two colliding particles must 
be comparable. If this is not the case, the denser cloud 
will likely cross the less dense one without feeling sizeable 
effects. Therefore, we assume that to have a collision the 
mean densities of the two colliding clouds should not differ 
by more than a factor of ten. Finally, depending on the 
mutual distance two kinds of collisions are possible: 

— The head-on collision, which happens if the distance 
between the clouds becomes shorter then their mean 
radius. In this case, the kinetic energies of the colliding 
particles are decreased to 20% of their original values, 
the directions of the colliding components of the ve- 
locity vectors are reversed, and the molecular density 
of both particles is increased by a factor of 100. 

— The off-centre collision, which happens if the distance 
between the clouds gets shorter than the sum of the 
two radii but larger than their mean radius, then their 
mean radius In this case, the molecular densities are 
increased by a factor of 10 and the kinetic energies are 
decreased to 80% of their initial value. 

— No interaction is supposed to occur in all other cir- 
cumstances. 

Clouds are allowed to experience collisions until they reach 
the molecular density of 10~^^g/cm ^ All assumptions 
are co nsistent with the results of the iMarinho fc Lepine 
(|200flf ) simulations. 

(6) The internal evolution of the clumps goes as fol- 
lows. Each "cold particle" represents a self-gravitating 
portion of the ISM in which two gas phases at differ- 
ent densities and temperatures coexist, i.e. the molecu- 
lar core of the cloud and the surrounding warmer mate- 
rial. We cannot determine the temperature and density of 
the two components, because physical processes at these 
scales are below the resolution power of our simulations. 
It is known that in cold collapsing clouds the " fragmenta- 
tion" process can break them into pieces when the density 
exceed the value set by the Jeans stability condition. As 
already anticipated, we cannot use the perfect gas equa- 
tion of state for the molecular part of the clouds, because 



E. Merlin & C. Chiosi; Multi-phase description of the ISM in early-type galaxies 



5 



different laws are expect to hold (i.e. lLarsonlll981l) and 
complicated processes such as turbulence and magnetic 
interactions may play very important roles in the molec- 
ular core. To partly get rid of this d ifficul ty, starting from 
the study of Springel fc HernquistI ( 20031 ) in cosmological 
context, we describe the cold cloud using mean values that 
refer to the "molecular part" and the "neutral" embed- 
ding warmer phase, and assume that in the cold particle, 
radiative cooling causes the increase of H2 phase at the 
expenses of the HI phase (the process is occurring at con- 
stant temperature for the two components) according to: 



5PH2 



5pHI ^Hl{PHI,UHl) 



Ot 01 UhI — 'UH2 

In practice, when a gas particle detaches itself from the 
surrounding hot ambient medium and becomes a " cloud" , 
it is thought to be made of molecular and neutral gas in 
some proportions (in our simulations, we assume that no 
H2 is present when the cloud is formed and it increases af- 
terwards). The radiative losses "cool" the neutral part and 
the "cooled" material is subtracted from it and accreted 
onto the molecular part, thus changing the densities of the 
two phases (at constant occupational volume). 

(7) Only cold particles are allowed to form stars. The 
transition from cold gas particle s to star particles is based, 
as explained in lLia et al.l (|200j) . on a probabilistic inter- 
pretation of the ISchmidd ( 19591 ) law. In practice, we con- 



sider that on a time scale of about 3 Myr (the lifetime of 
a typical high mass star) a fraction of the molecular gas, 
proportional to its density, is turned into stars according 
to 



dp^, 
~dt 



dpH2 

tdyn 



(3) 



Instead of creating more and more new particles of dif- 
ferent mass, which would make the computational cost too 
high, we prefer to interpret this formula as the probabil- 
ity that a cloud particle is completely turned into a star 
particle. The probability is 



P{SF) = 



PH2 



{PH2 + phi) 



1 — exp 



td' 



yn 



(4) 



where the first term takes into account that we consider 
only the m olecular density when t urning the whole particle 
into stars. Buonomo et al. I (|2000l) showed that, provided 
the total number of particles is sufficiently high (say more 
than ~ 2000), the right global behaviour is achieved (the 
statistical descriptio n equals the straight one). Following 
Lada fc Ladal (2003), we adopt c* = 0.05. As already men- 
tioned, close encounters between clouds enhance star for- 
mation, due to the growth of molecular densities after a 
collision. 

(8) Along with radiative cooling and energy dissipation 
by inelastic collisions, the thermal evolution of clouds is 
governed by the energy input by SN explosions. The SN 
energy feedback is perhaps the most difficult issue to ad- 
dress, because of its complex and poorly known physics. 



Furthermore, all physical processes concerning SN explo- 
sions and evolution of the SN remnants (SNR) are well 
below both our time and space resolution scales (namely, 
~ 1 Myr and ^ 1 kpc), so that we are forced to repre- 
sent them through "global" changes on large ISM regions. 
Here we adopt a simple description and try to minimize 
the number of free parameters. 

feedback energy from SN explosions is distributed 
among particles as follows. The "death" of star parti- 
cles and associated SN explosio ns are described with the 
statistical method prop osed bv lLia et all (|2002l ). see also 
iMerlin fc Chiosil (|2006h . In brief, a star particle (a SSP 
of a certain age) is sorted out through a Monte-Carlo 
method and its "death probability" is measured by the 
ratio between the amount of mass lost by stellar winds 
and SN explosions and total SSP mass. The probability 
goes from zero to one when the ratio goes from "zero to 
one" (from zero to a maximum value that ultimately de- 
pends for a given IMF on the mass interval within the 
SSP in which SN explosions are supposed to occur). In 
this way, we avoid instantaneous recycling and, at increas- 
ing number of star particles, we get closer and closer to a 
realistic description. When a star particle dies, its mass is 
completely turned back to the hot phase, and the energy 
feedback due to SN explosions within the SSP is spread 
over the particle itself and its neighbours, in the following 
way. First of all we compute the total gaseous mass rritot 
inside a sphere of radius equal to the smoothing length 
of the newly born gas particleH. rritot is thus given by the 
mass of the particle itself plus the mass of all cold clouds 
inside the sphere of radius equal to the smoothing length. 
The corresponding release of specific energy u [erg/g] is 



NsNEth 
rritot 



(5) 



where Nsn is the number of SN explosions within 

the dying parti cle, calculated according to the 

Greggio fc Renzinil (1983) rates, and Eth is the thermal 
energy produced by each explosion (see below). Each 
cloud is then given an amount of thermal energy cor- 
responding to its mass, decreased accordingly to the 
smoothing kernel of the newly formed gas particle, and 
the exceeding energy is given to the new gas particle 
and added to its original thermal energy. Furthermore, 
thermal energy is spread among the hot neighbour- 
ing partic les according to the n ormal SPH algorithm. 
Following iThornton et al.l (jl998[ ). we assume that the 
amount of thermal energy spread over the ISM by each 
SN explosion (which yields a total energy of 10^^ ergs 
per event) depends on the ambient gas conditions. Thus, 
we compute the mean values of density and metallicity 
rto and Z/Zq of the neighbouring particles, weighing 
each contribution on the smoothing kernel, and use the 
following relation to obtain the total amount of thermal 



^ Using the particle smoothing length instead of a fixed dis- 
tance causes feedback to be diflterently efficient in different en- 
vironments, and avoids the introduction of parameters. 
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enerRy to be spread into the surrounding medium 
Thoriaton et al.lll998l 



for details): 



Eth = 1.83 X Wn^'>-'\ZIZQ)- 
+ 1.23 X 1049n^"-24(Z/Zc3)-"-°« 
iilog[Z/ZQ\ > -2, and 



-0.24 



ergs, 



see 



(6) 



Eth = 5.53 X lO^^rio 
+ 1.78 X lO^^no 



-0.23 



ergs 



(7) 



otherwise. Finally, each cloud is given an amount of radial 
kinetic energy, still weighed on the kernel and computed 
in the same way as for the thermal energ y, using the em- 
pirica l relation Ekin = 1.69 x 10^" ergs of iThornton et al 
(|l998h . 



With this mechanism, clouds are allowed to "evapo- 
rate" back to the hot phase as soon as the temperature 
of their neutral part gets higher than a threshold value, 
that is naturally taken equal to the cloud formation tem- 
perature. In addition to this, we try to mimic the first 
phase of adiabatic expansion of the SNR and passage 
of the compressional front (see iMcKee fc: Ostrikeii Il977 , 
for details) by momentarily inhibiting radiative cooling in 
clouds which experience SN heating. By doing this we en- 
sure that sudden cooling of the particles cannot take place 
simply because soon after the explosion they happen to be 
immersed in a denser environment. 

(9) Fin a.lly, the chemica l evolution of the ISM is fol- 
lowed as in Lia et al. I (|2002h . It is worth noting, however, 
that cold particles are supposed not to experience chemi- 
cal enrichment, which is limited to the hot medium. Only 
when a cloud has been evaporated by SN energy feedback 
it can enrich its gas with metals. 

Fig. [2] shows a schematic summary of the different in- 
teractions among the baryonic components of a model 
galaxy with the new scheme. 

3. Testing the new recipe 

3.1. GalDyn, The Padova NB-TSPH code 

GalDyn is an evolutionary code for galactic systems, 
developed in Padova over the years. It is a fully la- 
grangia n NB-TSPH code, ba sed on the well known Tree- 
code b y [Barnes k. Hut ( 19861 ) and SPH algorithm bv lLucv 
(Il977h . For a detailed des c riptio n of the basic version o f 
the code s ee Carraro et al. J 1998 V Buonomo et al 
Lia et all (|2002l ): iMerhn fc Chiosil (|2006l ). For the pur- 
poses of this study, the code has been improved in several 
ways with respect to the previous versions. 

(1) First of all, an adaptive multiple time-step inte- 
gration scheme has been adopted. For each particle an 
individual time-step taken as the minimum between the 
acceleration time-scale, the velocity ti me-scale, and, for 
gas p articles, the Courant time-step (see lMerlin fc Chiosi 
20061 for details) is calculated using its current physical 
values. Within each time-step for the whole system, only 



Gatheiing of cold gas 
(radiative cooling) 




Fig. 2. The different interactions among the three bary- 
onic components of a galaxy with the multi-phase descrip- 
tion of the ISM. 



those particles with a sufficiently small individual time 
step are considered to be "active" and are moved on, us- 
ing a Drift-Kick-Drift Leapfrog integrator. This is similar 
to the implementation used by iSpringel et al., (,2001|) in 
Gadget 1. However, it is worth noting that the Drift- 
Kick-Drift algorithm may suffer non negligible problems 
for what concerns the conservation of the motion integrals 
(see Springe! 20051 ). Note that particles are not allowed to 
have a time-step lower than a minimum value (fixed equal 
to 10^^ of the Hubble time), and whenever required their 
smoothing length is changed to establish the right value. 

(2) The adopted smo othing kernel for SPH p a rticles , 
originally developed by iMonaghan fc Lattanziol ( 19851 ). 
has been revised following Thacker et al.l (|2000f ) and refer- 
ences therein. In order to avoid artificial clumping in high 
density regions, the first derivative of the kernel function 
is set constant and equal to l/{TTh^) for small radii. 

(3) Radiative cooling is not allowed to be active on 
a particle when its temperature falls below the threshold 
value given by 



= C 



Pg 



2/3 



2/3' 



Pa, (8) 



see 



Bonnor 19561 ) where 



Too 



is the mass of a baryonic 
particle, fc is the Boltzmann constant, /i and e are the 
mean molecular weight and the gas softening parameter, 
mil is the proton mass, and the C is a constant (in cgs 
units, C=0. 89553). Below this temperature the particle 
would become Jeans unstable (obviously, adiabatic cool- 
ing is not limited by this condition). This somehow im- 
poses an unavoidable limit to the mass resolution of the 
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particles because, if the threshold value is too high, cool- 
ing and cloud formation cannot be correctly modelled. In 
Fig.[3]the threshold temperature in Kelvin degrees is plot- 
ted versus the logarithm of density (in g/cm'^), for several 
values of the particle masses and softening lengths. The 
solid lines refers to the values adopted in the present sim- 
ulations, see Sect. 13.21 

(4) In order to correctly evaluate the impact of the 
inverse Compton cooling at high redshifts, we roughly es- 
timate the ionization fraction and thus the number density 
of free electrons as a function of temperature. The number 
density is 



He = UHI if T > 10'' K, 

rig — riHi X ea;p(— (10'* — T)) otherwise, (9) 

where the value of 10'* K roughly corresponds to the ion- 
ization temperature of Hydrogen in astrophysical environ- 
ments. 

(5) We describe the thermal conduction 
among SPH particles as o riginally developed by 
Monaghan Sz Lattanziol ( 1991 ). Thermal conduction 
is expected to be important in situations of shocks, and 
in particular after a SN event. iMcKee &: Ostriked (|l977l ) 
pointed out that thermal conduction could evaporate 
cold clouds and warm up the ISM inside a SNR (thus 
also changing the density inside the remnant). Although 
in our scheme clouds are evaporated "automatically" 
due to SN events, we expect thermal conduction to be 
important in the subsequent phases of evolution, when a 
large number of hot particles are present in a high density 
region. 

(6) To take into account UV flux and stellar winds 
from massive stars without too much computational ef- 



fort, we make use of the results of Buonomo et al. I (|2000h 



who found that the contribution to the energy feedback by 
these two sources of energy is comparable to that of SN. 
Since stellar winds and UV flux come the same stars that 
eventually will explode as Type II SN, we crudely lump 
together the three contributions and simply increase by a 
factor of 3 the energy injection by SN explosions. 

3.2. Formation of a galaxy in cosmological context 

We have tested the new prescription for star formation 
implemented in GalDyn by simulating the formation of 
a galaxy from quasi-cosmological initial conditions. 

The cosmological s cenario we adopt is consistent with 
the W-Map3 results ( Spergel et al. 2006h . namely ho ~ 
0.73, = 0.238, = 0, 0-8 0.74, n = 0.951, and initial 
baryon fraction ~ 17.6%. 

The ini tial conditions are ob tained using the Grafic2 
package bv lBertschingei (2001) to generate an initial per- 
turbed cubic grid of 64^=262,144 particles. In this grid we 
search a region containing an overdensity peak and sin- 
gle out another much smaller cubic volume around it, so 
that the peak is located near its geometrical centre. This 
small cubic region is initially made of 6'^ =2 16 particles 




Fig. 3. Lower threshold temperature for SPH particles of 
mass nip and softening length e, as a function of their 
density. Solid line: nip = too,p = 6 x 10^ Mq, e — eg — 1 



kpc. Dotted lines: from top to bottom, = mo,p and e 
— 0.1, 0.5, 2 times Eq. Dashed lines: from top to bottom, 
e — eo and nip — 2, 0.5 and 0.1 times too,p- 



(corresponding to a cubic volume of ~ 20'^ physical kpc^) 
and is then re-simulated with a finer resolution, which is 
imposed to give at the end a number of particles com- 
patible with our computing capability. Finally, in the re- 
simulated small cubic volume we isolate the sphere whose 
radius is half of the cube dimension and whose centre coin- 
cides with that of the cube (i.e., with the original position 
of the density peak in the large cubic grid). This final 
sphere contains ^ 14,000 particles (approximately half of 
dark matter and half of baryons). The sphere is inserted 
into the conformal Hubble flow with void boundary condi- 
tions. These latter may not be fully adequate to describe 
the real situation because the isolated volume is in real- 
ity immersed in a background exerting some dynamical 
effects on the matter contained inside the sphere under 
consideration. To cope with this, we add a small initial 
solid-body rotation to the sphere, and follow its expan- 
sion, detaching from the Hubble expansion , turn-around, 
and collapse (see e.g. lMerlin &: Chiosi 20061 . for details on 
the initial conditions). The sphere can be considered as 
a good approximation of an isolated initial proto-galaxy. 
Despite the small solid body rotation, the void boundary 
conditions may cause some probl ems, which anyway are 
not crucial for our purposes (see iMerhn & Chios] l2006l 
for more details). 

As already said, the way the cubic volume single out 
from original grid has been re-simulated is chosen in such 
a way that the total number of particles in the so-called 
proto-galaxy is compatible with our limited computing 
facilities. We are well aware that the rather small total 
number of particles may constitute a major drawback of 
the results we are going to present. The main motivation 
for this choice is that we want to test the new scheme 
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of star formation with an affordable comp uting cost and Table 1. Initial parameters for the model proto-galaxy 



time. 



However several studies (see e.g. Kawata 1999t 

Buonomo et ar200lj; 'Kawat J 2001al lbl: iKawata k Gibson 



2003; Scmclin & Combes 200 J have already shown that 
a total number of particles of the order of lO"* is enough 
to achieve a sufficiently accurate global description of the 
physical processes. Therefore, we are confident that the 
main body of the physical conclusions we are presenting 
are not too badly influenced by resolution limits. They 
have to be viewed as indications of the right trail to fol- 
low in the future. 

Wc also consider the action of the re-ionizing back- 
ground radiation (e.g. the UV field of the PopIII massive 
stars). To deal with effect in a very simple way, we crudely 
add a constant heating of 10""^"' ergs/s/atom from redshift 
z=20 to z=6. 

As the simulation starts well before re-ionization and 
the initial overdensity of the proto-galactic halo is small. 
Hydrogen can be considered as almost completely neutral 
at the initial stage. Anyway, the initial temperature of the 
gas is set equal to the maximum between a " cosmo l ogical " 
value, obtained as explained in Merlin fc Chiosi ([2006') , 
and the minimum temperature allowed by the Jeans limit 
discussed above. Indeed, the cosmological temperature of 
baryonic matter at z ~ 60 would be below 100 K, but the 
Jeans resolution in the present models imposes a thresh- 
old value of ~ 2000 K. As the Jeans threshold is well 
below the value at which ionization is expected to yield 
sizeable effects, we can safely use this value instead of the 
real cosmological limit. In any case, at the very initial 
stage we must consider the gaseous component as made 
only by "hot" particles, because as already explained, cold 
clouds are allowed to form only when belonging to collaps- 
ing flows. 



Ballesteros-Paredes et al.l (120071 ) estimated the follow- 



ing densities and temperatures for the two main compo- 
nents of a cold cloud: 1 < n < 100 cm"^ and 100 < T < 
5000 for the atomic gas; 100 < n < 1000 cm~^ (up to 
10^ in cores) and 10 < T < 100 K for the molecular gas. 
Therefore, we adopt Tcoid = 5000 K and pcoid = 10"^'' 
g/cm^ as threshold values. Such a low value for the den- 
sity threshold is motivated by the fact that this is the 
global SPH density of a particle, while in order to create 
cold clouds we are only interested in the densest part of it. 
We then choose acoid — 1000 K. Finally, to describe the 
properties of of our star particles (energy feedback and 
chemical enrichment) we adopt the Salp eter (1955) IMF. 

Our model proto-galaxy has an initial radius of 19 kpc, 
and a total (Dark and Baryonic) mass of 2.5 x lO^^M©, 
with an overdensity of ~ 0.09 at the initial redshift z = 60. 
The number of particles for Baryonic and Dark Matter is 
~ 7000, the same for each component. For baryons the 
mass of a particle is ~ 6 x IO^Mq, and the softening length 
is 1 kpc; the corresponding values for the Dark Matter 
component are ~ 3 x lO^M© and 2 kpc, because of the 
much higher total mass. All the relevant quantities in use 
are summarized in Table [TJ 



Initial redshift 


60 


Initial mean overdensity 


0.09 


Initial radius [kpc] 


19 


Total mass [Mq] 


2.5 X 10'' 


Number of baryonic particles 


7000 


Mass of a baryonic particle [Mq] 


6 X 10" 


Softening length for a baryonic particle [kpc] 


1 


Number of Dark particles 


7000 


Mass of a Dark particle [Mq] 


3 X 10' 


Softening length for a Dark particle [kpc] 


2 



-0.1 -0.05 0.05 0.1 -0.1 -0.05 0.05 0.1 



-0.1 -0.05 0.05 0.1 -0.1 -0.05 0.05 0.1 

Fig. 4. From top to bottom, from left to right: position 
(physical Mpc) of Dark Matter particles projected on the 
XY plane at redshift z = 5.9, 3.4, 1.7, 1.0. 
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Fig. 5. Same as above, but for hot gas particles. 
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Fig. 6. Same as above, but for cold gas particles. 



-0.1 -0.05 0.05 0.1 



-0.1 -0.05 0.05 0.1 



Fig. 7. Same as above, but for star particles. 

3.3. Formation and evolution of a galaxy: results and 
comments 

First of all, we note that owing to the present procedure 
to isolate the proto-galaxy from the cosmological context, 
the void boundary approximation may influence the evo- 
lutionary time scales and the detailed results in turn. 
In particular, it may (slightly) change the relationship 
between rest-frame ev olutionary time and redshift (see 
Merlin fc Chiosi 20061 ). Despite this cautionary remark, 
we are rather confident that the basic features and prop- 
erties of the model should not depend very much from this 
drawback of the present formulation. 

Galaxy formation. In Figs. |4]through[7]we show the 
spatial evolution of the system, in its four components: 
Dark Matter, hot gas, cold clouds and stars. All boxes are 
100 kpc per side. The number of cold clouds (whose di- 
mensions have been magnified for clarity) is always quite 
small, as they soon evolve into stars or are evaporated back 




3 3.5 

logLlJ [MyrJ 



Fig. 8. Fractional mass in the three baryonic components 
versus time. Dashed: hot gas; dotted: cold clouds; solid: 
stars. 



10 5 3 2 




Fig. 9. Star formation rate in solar masses per year. 

to the hot phase; this is clearly visible in Fig. [8] showing 
the fractionary mass of the baryonic components as a func- 
tion of the cosmic time. This situation, which differs from 
the one in a disk galaxy like the Milky Way, where most 
of the gaseous mass is stored in the cold phase, can be 
explained as follows. First, there is a technical reason: as 
already noted, at the beginning of the simulation, in order 
to apply the SPH formalism and get the right values for 
densities and temperatures, we consider the gaseous com- 
ponent to be made only by hot particles, even if their tem- 
perature is below the threshold value for cloud formation. 
This effect should become less relevant as the galaxy ages. 
Second, the mass resolution is small, so that statically the 
number of cold clouds going into stars per time step is also 
small (and fluctuating). The third reason is more physi- 
cal: in our model galaxy (simulating what happens in an 
elliptical) a large fraction of gas is turned into stars so 
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that little gas is left after a few Gyr of star formation ac- 
tivity (see below). Most of the cold gas sinks very early 
toward the centre of the proto-galaxy, whereby undergo- 
ing frequent collisions is turned into stars at a high rate. 
Dying stars soon provide SN feedback to the surrounding 
medium, thus inhibiting the formation of a great number 
of new cold clouds. Conversely, in situations with low star 
formation rates (such as in evolved disk galaxies) the cold 
gas is expected to distribute almost uniformly across the 
disk and likely to co-rotate with the rest of the galaxy; 
the number of cloud-cloud collisions is of course lower and 
lots of cold gas survive to star formation. See below for 
some additional comments on this issue. 

The clumps of collapsed matter merge at very early 
times, and at z ~ 3.5 the main body of the final object is 
in place, accreting new material from the outskirts. In the 
context of current terminology in literature, we prefer to 
classify this situation as "monolithic" rather than "hier- 
archical" . In any case it would be a " very early hierarchi- 
cal" assembly of a galax}|f|. At the end of the simulation 
(z ~ 0) , the galaxy has a nearly spherical shape with axial 
ratios h/a = 1.02 and c/a — 1.12. The total stellar mass 
is ~ 3 X 10^" Af© (75% of the initial baryonic mass), with 
an half mass radius i?i/2 = 4.6 kpc. Almost no gas is left 
inside the virial radius of the object. 

In Fig. [9] the rate of star formation (in solar masses 
per year) is shown. The first stars begin to form around 
z ~ 12. The bulk of stars (80% of the total) is in place at 
z ~ 1.6. The galaxy has already acquired its final features, 
and from now on it evolves almost passively. The last, mi- 
nor burst of star formation takes place at z ~ 0.7. These 
result is in good agreement with the obser vations for ob- 
jects with similar total stellar mass (see e.g. 
I2OO5I: IClemens et all l2006t ICimatti et aJ 



Thomas et al, 
2006|). In the 



nowadays well-established downsizing scenario, smaller 
objects assemble their stellar content over longer time 
scales than massive galaxi es - at least in high de nsity envi- 
ronments. Among others, Jimenez et al. ( 20051 ) studied a 
large sample of ~ 10^ galaxies from the SDSS, and found 
a tight relation between the total stellar mass of an ob- 
ject and epoch of its assembly (see e.g. their Fig. 4). A 
system with 10^" < M^,/Mq < 10^^ is expected to form 
most of its stellar content about ~ 5 — 6 Gyr ago, whilst 
more massive objects are expected to exhaust their activ- 
ity much earlier (the situation being more intrigued for 
low mass objects, which have different histories depend- 
ing on their initial overdensity). Our results are in very 
good agreement with this picture, and no particular fine 
tuning of parameters was required to obtain them. Indeed, 
all the parameters are conservative and taken from previ- 
ous studies of both observational and theoretical nature. 



^ We name monolithic a formation process in which all the 
relevant events which give a galaxy of any size its final shape 
and features, e.g. mergers and dynamical evolution, have taken 
place in the remote past, say before z ~ 2, in contrast to a 
hierarchical view in which bigger objects are assembled from 
smaller objects over the whole Hubble time. 



Perhaps, a more prolonged period of activity could have 
been expected for our object, given its mass and its rather 
low initial overdensity, but of course we cannot too strictly 
compare the results for a single model galaxy with the sta- 
tistically averaged trends that are found observationally. 

In order to achieve some insights on this point, as well 
as on the improvements due to the new prescription for 
star formation, we compared the present results with those 
obtained with the clas sical SPH treatment of the ISM (see 
Merlin fc Chiosi 20061 for details). The two star forma- 
tion histories are shown in Fig. [TOl where the solid lines 
are polynomial best fits of the numerical results, sketched 
in order to give an idea of the trends of the star formation 
histories in the two models apart from the very fluctuat- 
ing profiles, which are due partly to numerical artifacts 
and partly to the statistical method we have adopted to 
deal with star formation and energy injection by dying 
stars. The SPH simulation was stopped after 7 Gyr of 
evolution due to limited computational resources. With 
the new prescription, the formation of the stellar mass 
is more gradual and lasts longer; the peak of activity is 
lowered from ~ 17 to ~ 12 solar masses per year, and 
shifted from z ~ 4.3 to z ~ 2.9. The total stellar mass as- 
sembled in the two models is, however, almost equal. The 
different behaviour obtained with the two star formation 
prescriptions can be understood first of all by looking at 
the different time scales involved in the formation of new 
stars. In the classical SPH treatment, all gas particles can 
in principle be turned into stars, and therefore a large 
reservoir of fuel is ready to be used right from the begin- 
ning, giving rise to a strong and sudden burst of activity. 
In contrast, with the new scheme, the gas has first to gen- 
erate cold clouds, than these must develop a molecular 
core, and finally the core must reach a high density be- 
fore a star particle can be born (only the molecular part 
of the cold particle is considered in the statistical algo- 
rithm adopted to form stars). Another difference between 
the two models is the amount of energy injected into the 
gaseous component. While in the classical SPH treatment 
only the SN feedback was taken into account, in the new 
recipe we also considered other sources such as UV flux 
from young stars, stellar winds from massive stars, and the 
kinetic kick to cold clouds by SN explosions. As widely dis- 
cussed below, this has certainly sizeable effects on the star 
formation histories of the two models. All this end up in 
a longer and shallower period of activity. Of course, being 
the law of star formation at work always the same but for 
the amount of material to its disposal, the change in the 
global time variation is not dramatic, though sizeable. 

Tails of star formation. The model shows a very 
long tail of stellar activity in the very central regions which 
might be taken as a fl aw to be eliminated. However, some 



recen t studies (see e.g. iBressan et al.ll2006t IClemens et al 



20061 ) argue for the presence of a small tail of ongoing 
star formation up to low redshifts, at least in clusters, for 
a conspicuous fraction of early type galaxies. With the 
present model, whether the tail corresponds to that sug- 
gested by the observations cannot be said because of the 
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Fig. 10. Star formation histories. Dotted: the new multi- 
phase description of the ISM; dashed: the standard SPH 
model. The solid lines are best fitting polynomial func- 
tions, see text for details. 



[kpc] 



Fig. 11. Surface density profile of star particles. The solid 
line is a fitting Sersic profile (m = 3). 

the mass resolution. Indeed, if the number of gas particle 
turned into stars during the tail activity is very small (of 
the order of one to two particles per time step) the mass 
of newly created stars is always significant because of the 
mass resolution. We expect that increasing the mass reso- 
lution (initial total number of particles) should clarify the 
issue and lead to quantitative estimates of the star for- 
mation efficiency in the tail activity (if any). On the other 
hand, a mechanism able to completely stop star formation 
in most massive early type galaxies is still missing. Many 
hints suggest that the presence of an AGN could be the 
solution. Work is in progress to address the problem. 

Finally, we would like to note that the mass in stars of 
ou r model is almost equal t o the limit of 3 x lO^^MfT^ f ound 
bv lKauffmann et al. and I Jimenez et al.l (|2005l ). be- 



Fig. 12. Mean chemical composition versus time. Solid: 
stars; dashed: hot gas; dotted: stars with the old descrip- 
tion of the ISM. 
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Fig. 13. Metallicity gradient for star particles. 



tween "red and massive" and "blue and small" objects in 
the SDSS. In this sense, we could expect to obtain an ob- 
ject with intermediate features between the two galactic 
populations: the star formation history of our model goes 
in the right direction. 

Surface mass density profile. In Fig. [Tl] we show 
the surface de nsity profile o f the final stellar component 
together with a Sersic ( 1968h profile with m = 3. The pro- 
file fits quite well the Sersic relation, especially if we con- 
sider that our profile is actually a projected mass profile, 
which should be photometrically translated into a lumi- 
nosity profile in order to obtain a better comparison with 
the observational profiles. 

Chemical enrichment. Fig[T2] presents the metal- 
licities of stars and hot gas. The final value for the star 
metallicity, at z ~ 0, is solar within a few percent, i n strik - 
ing agreement with the results by iGallazzi et al. I (|2005l) . 
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By contrast, in the "old" model the metallicity is already 
super-solar at z ~ 4, and at z ~ 0.7 is 1.3 Zq, which 
is expected for objects ten times as massive (dotted line 
in Fig. [T^ . Likely, this is partly due to the same mecha- 
nism that decreased the star formation rate, and partly to 
the procedure adopted to spread metals into the interstel- 
lar medium. As already explained, metals are indeed dis- 
tributed only among hot gas particles, leaving cold clouds 
with their own original metallicity until they are heated 
up back to the hot phase by SN feedback, thus delaying 
the chemical enrichment of the material out of which stars 
are formed. 

Metallicity gradients. Fig. [13] shows the radial 
metallicity gradient of the stellar component. It has been 
long known t hat elliptical galaxies are optically redder 
in their cores ( de Vaucouleursl[l96l[ ): this has been inter- 
preted in terms of higher metallicity due to prolonged star 
formation in the central regions. This is well achieved in 
our model. 

Galactic winds. Finally, conspicuous galactic winds 
are found to occur. The upper panel of Fig. [Ml shows the 
mean radial velocity of the gas particles at z ~ 0, as a 
function of their radial distance, compared to the escape 
velocity at the same distance. For the sake of comparison, 
the lower panel shows the same plot for an identical model 
where energy feedback had been drastically decreased (al- 
most switched off, see below) : the number of gas particles 
are fewer (one third of total gaseous mass), and almost all 
of them are confined within the galaxy potential well. To 
avoid possible misunderstanding of the physical meaning 
of Fig.[T3|it is worth commenting on the particles that are 
found at very large distances from the mother galaxy. This 
is a result of the void boundary conditions and the com- 
putational scheme in use. Particles that acquire velocities 
larger than the escape velocity and are free to travel with- 
out feeling other effects can indeed reach very large dis- 
tances (for instance a particle escaping at about 400 km/s 
and keeping this velocity constant for about 10 Gyr travels 
up to ~ 4 Mpc). Second the numerical code keeps track 
also of these particles even if they do not longer belong to 
the galaxy. Fig. [T3| simply shows that large amounts of gas 
can be considered unbound to the galaxy, because of the 
energy feedback, and transferred to the external environ- 
ment as a galactic wind. This has slightly super-solar mean 
chemical composition and will affect the chemical com- 
position of the medium in which galaxies are immersed, 
the well known issue of the chemical composition of the 
intra-cluster medium (see Chiosillioool : Moretti et al.ll2003t 



Portinari et al. 2004L for recent discussions of this sub 



ject). The results for the model with drastically decreased 
energy input also deserves some comments. In this mod- 
els the energy input is not completely switched off because 
otherwise all gas particles around would have been turned 
into stars and no gas would have been left over. The energy 
decrease is modulated in such a way that the energy input 
brings gas particles to the temperature of about 50,000 K, 
i.e. half of the hydrogen first ionization potential (13.6 eV 
per atom) and well above the temperature at which the 




Fig. 14. Upper panel: mean radial velocity of gas particles 
versus radial distance, compared to the escape velocity 
at the same distance (dotted line). Lower panel: same as 
above, for a model with SN feedback nearly switched off. 



astrophysical low-density plasma starts ionizing and be- 
low which it begins to form cold clouds. Compared to the 
case with full energy input, in which the temperature of 
the gas particles can reach some 10'' — 10® K, there is a 
difference of at least a factor of 100. The kinetic energy 
acquired by the gas particles scales with the square root 
of the temperature and therefore the cloud velocities are 
now about a factor of 10 smaller than in the full energy 
input cases. Over a time interval of 10 Gyr these particles 
can travel about 1 Mpc and some of them can even acquire 
sufRcient energy to exceed the escape velocity, consistently 
with the results shown in Fig. [TH 

Star formation and energy feedback. The rela- 
tionship between star formation rate and energy feed back 
deserves particular attention and some more detailed in- 
vestigation. In general, one expects that the higher the 
energy given to the interstellar medium by SN-explosions, 
UV radiation, stellar winds, the lower is the efficiency of 
star formation. This holds in particular with our scheme 
in which only cold gas is allowed to form stars. Our ref- 
erence model (the one we have presented) contains four 
sources of energy feedback: SN-explosions, UV radiation, 
stellar winds, and a contribution of radial kinetic energy 
imparted to gas clouds surrounding a dying star particle 
suffering SN-explosions (see Sec. 13.11 point 6). This sim- 
ple scheme could be too a crude description of reality. In 
addition to that, our detailed history of star formation 
proceeds in a number of short bursts of activity, which 
could be ascribed partly to the scheme we have adopted 
to deal with energy feedback and partly to numerical reso- 
lution (this latter point will be examined below) . In order 
to investigate the effect of energy feedback on the star for- 
mation rate, we run four simulations with the same initial 
conditions and total mass of the galaxy model we have 
being discussing in so far. These four models are, how- 
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ever, calculated with 3500 particles of baryonic mass and 
3500 particles of dark matter, i.e. half of the number for 
each component contained in the original model. The rea- 
son of this choice will become clear when discussing the 
effect of the numerical resolution. The difference among 
the four models is only in the energy feedback: (1) Model 
A is the same as the original model but for the fewer 
number of particles. It contains all sources of energy feed- 
back. (2) Model B contains only the SN-explosions, UV 
radiation and stellar winds (no contribution by radial ki- 
netic energy to nearby gas particle). (3) Model C contains 
only the SN-explosion and contribution of kinetic energy 
to nearby gas particles (no UV and stellar winds inputs 
are considered). (4) Finally, model D has only the contri- 
bution from SN-explosions, while all the rest is quenched 
off. Fig. [15] compares the results for star formation (as 
in Fig. [TOl in this case too polynomial best fits of the 
star formation rate are superposed to the fluctuating star 
formation rate of the numerical model to highlight the 
general trends). Fig. [TBI compares the results for chemical 
enrichment in the four models. As expected, the effects 
are indeed remarkable. Including more and more sources 
of energy feedback tends to lower the mean star formation 
activity and to quench it much earlier as compared to the 
cases in which they are partially included. Consequently, 
chemical enrichment also stops and saturates much earlier. 

The separated contribution by the different sources 
of energy feedback is best illustrated in Fig. [17] show- 
ing the cumulative building up of the stellar mass for the 
four models. Models A and B contain all thermal energy 
sources (SN, UV, Winds) and differ for the kinetic kick. 
The same applies to the pair C and D, which however 
contains only the SN input for the thermal part. Little or 
no difference is caused by the kinetic energy kick (com- 
pare A with B and C with D). The largest difference in 
the cumulative stellar mass is due to the differences in the 
thermal energy input: star formation efficiency and hence 
total stellar mass decrease at increasing feedback. 

The star formation history always shows a high num- 
ber of burst episodes, whose amplitude seems to in- 
crease at decreasing the number of energy feedback terms. 
However, this can be partly ascribed to the lower number 
of particles used in these simulations, as discussed next. 

Numerical resolution and physical results. What 
is the effect of the number of particles (mass resolution) 
used in our models on the physical consistency of the 
results? This is an important issue to address consid- 
ering that multi-phase NB-TSPH simulations of galaxy 
formation and evolution with significantly higher num- 
bers of part ic les ca n be found in literature: for instance 
from 20,000 to 250,000 particles). 



Booth et all EoO 



Scaiinapieco et al 



(l2006l from 18,000 to 80,000 particles), 
and iHarfst et al.l (12006 . 35,000 particles). These studies 



present also higher number of cases, whereas we present 
only one case. Therefore, compared to those studies our 
resolution and number of models are somewhat limited. 
As far as the number of models is concerned we put it 
aside because a statistical analysis of the variety of pos- 




Fig. 15. Star formation histories. Top left panel: model A; 
top right: model B; bottom left: model C; bottom right: 
model D. See text for details. 
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Fig. 16. Chemical enrichment histories. Top left panel: 
model A; top right: model B; bottom left: model C; bot- 
tom right: model D. See text for details. 



sible results (model galaxies) is beyond the aims of the 
present study. The crucial question to address is whether 
our resolution is sufficient to yield physically consistent re- 
sults. We have already mentioned that a number of studies 
have been made based on a rather small number o f parti- 
cles comparable to ours, e.g. iBuonomo et~al] (l200d 10000 
of dar k Matter and 10000 of Baryonic Matter) . iBuonorno 
(I1999L who lowered the p articles of Baryonic Matter down 
to a few thousands) and Semelin fc Comb^ ( 20021 who 
used 50000 particles of Baryonic Matter and 10000 of Dark 
Matter). All these models adopted "isolated initial con- 
ditions". Other studies adopt cosmological initial condi- 
tio ns and numbe r s of p art icles comparable to our s (see 
e.g. lKawata|[l999l l2001allbl : iKawata fc Gibsonll2003l ). 
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Fig. 18. Comparison of the star formation histories (top 
panel) for the original model with 14,000 particles (about 
7000 of Dark Matter and 7000 of Baryonic Matters) and 
full energy input indicated by the long dashed line and 
model A with the same energy input but 7000 particles in 
total (about 3500 of Dark Matter and 3500 of Baryonic 
Matter) indicated by the dotted line. The polynomial fits 
to star formation rate are shown for both cases to guide the 
eye in the gross evaluation of the star formation intensity 
and its temporal evolution. The bottom panel show the 
same but for the chemical enrichment history 

Star Formation. It is worth comparing the star for- 
mation and chemical enrichment histories of our fiducial 
model with those of case A. The models have exactly the 
same initial conditions, total mass, and energy input, but 
they differ in the number of particles (14000 versus 7000 
in total). This is shown in the top and bottom panels of 
Fig. [181 In the top panel of the star formation rate we 



also plot the polynomial best fits to guide the visual in- 
spection of the results. The overall trend is the same, even 
if the star formation rate of model A is more fluctuating 
and slightly less intense than that of the original model. 
The difference at the activity peak is ^ 23%, which is 
immediately mirrored by the chemical enrichment history 
(bottom panel) where the maximum metallicity is ^ 15% 
less than in the original model. These effects are clearly 
caused by the number of particles. As far as the fluctua- 
tions in the star formation rate are concerned, there are 
at least three causes: (i) the stochastic treatment of star 
formation rate and energy input by SN explosions, stellar 
winds and UV radiation; (ii) the small number of par- 
ticles (in particular those prone to star formation) and 
(iii) the age sampling of the numerical results, which are 
displayed averaging the numerical results over time in- 
tervals of 50 Myr. Certainly, by increasing the number 
of particles and in particular those of cold molecular gas 
(see below), they should smoothen out, if not disappear. 
Therefore we dot attribute them any particular physical 
meaning. What matters here if the global trend repre- 
sented by the analytical fits. One may guess that increas- 
ing the number of particles to about 50000 would be the 
ideal situation. The parallel version of the same code is 
in the final testing phase (at CINECA Super-Computing 
Centre) , so that simulations with much higher numbers of 
particles are next in our plans. 

Number of molecular clouds. At each time step, the 
number of cold gas particles in the fiducial model is of 
order of one hundred, which leads to a strong Poisson 
noise for star formation rate, which is enhanced by cloud- 
cloud collisions. Does it reflect onto the structural proper- 
ties, e.g. the Sersic-like mass density profile of our model 
galaxy? In Section 3.3 we have provided a plausible expla- 
nation why in the physical situation of our models aimed 
at simulating elliptical galaxies the number of cold clouds 
remains small as compared to that expected in a disc 
galaxy. It will not be repeated here. In any case, the num- 
ber of cold gas particles is self-consistent with the physical 
input of the models. Increasing it further is not possible. 
Whether increasing the number of particles, we would end 
up with a disc galaxy instead of an elliptical one is the 
subject of the next paragraph. 

Disc or elliptical galaxies? The models we have calcu- 
lated (the fiducial one in particular) strictly resemble an 
elliptical galaxy (as far as we can tell looking at the struc- 
tural parameters, surface mass density profile and axial 
ratios). Can a simulation, calculated with a larger number 
of particles but the same initial conditions and physical in- 
put, end up as disc galaxy? We do not think so. Very little 
angular momentum is initially given to the proto-galaxy, 
and there are no external torques which could induce size- 
able rotational effects. Therefore, the formation of an al- 
most spherical structure is a straight consequence. It has 
a reasonable density profile and a stable relaxed config- 
uration which resembles real spherical/elliptical galaxies. 
Likely perturbations with little initial rotation and weak 
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external torques will end up with a nearly spherical shape 
and Sersic profile of the mass distribution. 

Finally, similar results h ave been presented by (see 
e.g. Kawata 1999, 2001a, b; Ka wata fc G ibson 2003) using 
similar initial conditions, number of particles, and input 
physics, but for the multi-phase description of interstellar 
medium and star formation efficiency. 

4. Discussion and conclusions 

The new description of the ISM allows us to get a deeper 
insight on the time scales and processes of star formation, 
even at the relatively low resolution of the present model. 
Good agreement with important observational informa- 
tion is achieved, especially for what concerns the mass as- 
sembly time scales and chemical evolution of the system. 
The agreement can be further improved by finely tuning 
some of the model parameters. 

More and more observational evidences hint the ex- 
istence of a population of big, red, and isolated galaxies 
already in place up to z ~ 6 (see e.g. Mobasher et al.llioOSt 
Cimatti et al. 19991 20041 ). Moreover, recent large sur- 



veys like DEEP2, SDSS etc. clearly demonstrate that not 
only the stellar populations of massive galaxies are older 
than the one s belonging to smaller systems {downsizing, 
iBundv et al., 2005jf| but that massive galaxies assembled 
their mass at higher r edshift than the les s massive ones 
(top- down formation, Bundv et al.l l2006t ICimatti et al 



It is soon clearlI^^^^downs^z^ng and top- 



down views of star formation and mass assembly are in 
conflict with the "classical" predictions of the hierarchi- 
cal theories of structure formation, whereas they can be 
easily framed into the monolithic (as explained above) 
view of galaxy formation. Despite the many attempts 
that are being made to reconcile the hierarchical sce- 
nario with the new compelling evidences of "very early 
and nearly mo nolithic " assembly of early type galaxies 

- dry mergers (Bell et al. 2004) and hierarchical down- 
sizing ( De Lucia et al.ll2006h are the most popular ones 

- the observational data leave little room for viable so- 
lutions in the context of the hierarchical scheme. In par- 
ticular, we note that dry mergers cannot have played a 
fundamental role in the formation of early- type galaxies, 
given the weak dependence of their formation time-scales 
on environmental conditions and weak cosmic evolution of 



the mass funct i on of galaxies ( Bundv et al.ll2006h . Finally, 
Clemens et all (|2006[ ) show how a burst of monolithic star 



formation in high mass galaxies and a slower, less efficient, 
ongoing star formation activity in smaller objects are nec- 
essary conditions to achieve the observed trends of metal- 
licity and a- enhancement. Remarkably, down-sizing and 
top-down are natural products of the models presented by 



The c oncept of downsizing star formation was first sug- 
gested by iMatteuccil ([l994) from v ariations of the [Fe/M g] ra- 



tio with the galaxy luminosity, by iBressan et al.1 (ll996D from 
th e analysis of the lin e absorption indices and UV excess, and 
by lCowie et all (|l996l ). 



Chiosi fc Carraro ( 20021 ) (see they star formation histories 
as a function of the total mass and initial density) . 

A monolithic-like star formation as the one we 
present here (as wel l as in previous studies, see e.g. 
Merlin fc Chiosill2006f) . with more or less prolonged du- 
rations of star formation activity (depending on the total 
mass and initial density of the proto-galaxy, and perhaps 
on the environment in which it is being formed) seems 
to be able to match most of the characteristics of old, 
massive spheroids, as well as of younger and bluer ob- 
jects. We plan to perform an extended study on the de- 
pendence of the star formation history on the mass and 
density of the proto-galaxies, using the parallel version of 
GalDyn in which the new prescription is being imple- 
mented. The longer period of activity found in the new 
models agrees with the trends suggested by the observa- 
tional data. Fine tuning of the model parameters is the 
first step to undertake to reproduce the pattern of ob- 
servational constrains. The goal is to account for key di- 
agnostic planes suc h as the classical Fundame ntal Plane 
for local ellipticals ( Diorgovski fc Davis 19871). the mass 
quenching relation found bvlBundv et al. I (l200fih . and the 
chemical trends found by iGallazzi et all (|200ft . 

Given this picture, it is no longer necessary to consider 
mergers between proto-galaxies (or disks) as the main way 
in which massive ellipticals are formed. It is most likely, 
indeed, that this scenario has to be definitely discarded. 
If so, a general scenario reconciling our understanding of 
matter aggregation in the Universe on large scales with 
the compelling evidences about the formation of galaxies is 
missing. Likely, Nature follows the hierarchical mode when 
aggregating matter on the scale of groups and clusters, 
and the monolithic-like mode when aggregating (baryonic) 
matter on the scale of individual galaxies. On the other 
hand galaxy mergers cannot be completely ruled out, sim- 
ply because we have direct observational evidence of this 
phenomenon. They are beautiful, spectacular events, but 
not the dominant mechanism by which galaxies (the early- 
type in particular) are assembled and their main features 
imprinted. In the forest of galaxy formation theories, ex 
pluribus unum. 
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